https://rspatial.org/raster/rs/2-exploration.html
unzip('data/rsdata.zip', exdir='data')
## Warning in unzip("data/rsdata.zip", exdir = "data"): error 1 in extracting from
## zip file
“Create RasterLayer objects for single Landsat layers (bands)”
library(raster)
## Loading required package: sp
# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')
see the spatial resolution, extent, number of layers, coordinate reference system, etc
b2
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B2.tif
## names : LC08_044034_20170614_B2
## values : 0.0748399, 0.7177562 (min, max)
b3
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B3.tif
## names : LC08_044034_20170614_B3
## values : 0.04259216, 0.6924697 (min, max)
b4
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B4.tif
## names : LC08_044034_20170614_B4
## values : 0.02084067, 0.7861769 (min, max)
b5
## class : RasterLayer
## dimensions : 1245, 1497, 1863765 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : LC08_044034_20170614_B5.tif
## names : LC08_044034_20170614_B5
## values : 0.0008457669, 1.012432 (min, max)
# coordinate reference system (CRS)
crs(b2)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b3)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b4)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b5)
## CRS arguments:
## +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497 1
# spatial resolution
res(b2)
## [1] 30 30
# Number of bands
nlayers(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)
## [1] TRUE
s <- stack(b5, b4, b3)
s
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3
## min values : 0.0008457669, 0.0208406653, 0.0425921641
## max values : 1.0124315, 0.7861769, 0.6924697
plot(s)
plot(b5)
# create list of raster objects
filenames <- paste0("data/rs/LC08_044034_20170614_B", 1:11, ".tif")
filenames
## [1] "data/rs/LC08_044034_20170614_B1.tif"
## [2] "data/rs/LC08_044034_20170614_B2.tif"
## [3] "data/rs/LC08_044034_20170614_B3.tif"
## [4] "data/rs/LC08_044034_20170614_B4.tif"
## [5] "data/rs/LC08_044034_20170614_B5.tif"
## [6] "data/rs/LC08_044034_20170614_B6.tif"
## [7] "data/rs/LC08_044034_20170614_B7.tif"
## [8] "data/rs/LC08_044034_20170614_B8.tif"
## [9] "data/rs/LC08_044034_20170614_B9.tif"
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"
“The layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
We won’t use the last four layers and you will see how to remove those in following sections.”
landsat <- stack(filenames)
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 11 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values : 9.641791e-02, 7.483990e-02, 4.259216e-02, 2.084067e-02, 8.457669e-04, -7.872183e-03, -5.052945e-03, 3.931751e-02, -4.337332e-04, 2.897978e+02, 2.885000e+02
## max values : 0.73462820, 0.71775615, 0.69246972, 0.78617686, 1.01243150, 1.04320455, 1.11793602, 0.82673049, 0.03547901, 322.43139648, 317.99530029
| # Maps ## Grey-scale Single band and composite maps plot individual layers in a 2x2 grid legend values between 0-1 “different surface features reflect the incident solar radiation differently. Each layer represent how much incident solar radiation is reflected for a particular wavelength range. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.” |
r par(mfrow = c(2, 2)) plot(b2, main = "Blue", col = gray(0:100 / 100)) plot(b3, main = "Green", col = gray(0:100 / 100)) plot(b4, main = "Red", col = gray(0:100 / 100)) plot(b5, main = "NIR", col = gray(0:100 / 100)) |
| ## True / natural colour image “(vegetation in green, water blue etc), we need bands in the red, green and blue regions. For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used.” |
r landsatRGB <- stack(b4, b3, b2) plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite") |
| ## Compare True and False Colour Images |
| ```r par(mfrow = c(1, 2)) plotRGB(landsatRGB, axes = TRUE, stretch = “lin”, main = “Landsat True Colour Composite”) |
| landsatFCC <- stack(b5, b4, b3) plotRGB(landsatFCC, axes = TRUE, stretch = “lin”, main = “Landsat False Colour Composite”) ``` |
help(plotRGB)
# select first 3 bands
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 11 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values : 9.641791e-02, 7.483990e-02, 4.259216e-02, 2.084067e-02, 8.457669e-04, -7.872183e-03, -5.052945e-03, 3.931751e-02, -4.337332e-04, 2.897978e+02, 2.885000e+02
## max values : 0.73462820, 0.71775615, 0.69246972, 0.78617686, 1.01243150, 1.04320455, 1.11793602, 0.82673049, 0.03547901, 322.43139648, 317.99530029
landsatsub1 <- subset(landsat, 4:2)
landsatsub1
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2
## min values : 0.02084067, 0.04259216, 0.07483990
## max values : 0.7861769, 0.6924697, 0.7177562
# same as
landsatsub2 <- landsat[[4:2]]
landsatsub2
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2
## min values : 0.02084067, 0.04259216, 0.07483990
## max values : 0.7861769, 0.6924697, 0.7177562
# Number of bands in the original and new data
nlayers(landsat)
## [1] 11
nlayers(landsatsub1)
## [1] 3
nlayers(landsatsub2)
## [1] 3
landsat <- subset(landsat, 1:7)
landsat
## class : RasterStack
## dimensions : 1245, 1497, 1863765, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7
## min values : 0.0964179114, 0.0748399049, 0.0425921641, 0.0208406653, 0.0008457669, -0.0078721829, -0.0050529451
## max values : 0.7346282, 0.7177562, 0.6924697, 0.7861769, 1.0124315, 1.0432045, 1.1179360
names(landsat)
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)
## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"
“used to limit analysis to a geographic subset of the image”
extent(landsat)
## class : Extent
## xmin : 594090
## xmax : 639000
## ymin : 4190190
## ymax : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
e
## class : Extent
## xmin : 624387
## xmax : 635752
## ymin : 4200047
## ymax : 4210939
# crop landsat by the extent
landsatcrop <- crop(landsat, e)
landsatcrop
## class : RasterBrick
## dimensions : 363, 379, 137577, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 624390, 635760, 4200060, 4210950 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## source : memory
## names : ultra.blue, blue, green, red, NIR, SWIR1, SWIR2
## min values : 0.101796150, 0.075989284, 0.044045158, 0.030556191, 0.007243267, 0.001865030, 0.003144530
## max values : 0.4548080, 0.4746728, 0.5034941, 0.5592065, 0.7013395, 0.9038041, 0.9750224
par(mfrow = c(1, 2))
landsatsub3 <- subset(landsat, 4:2)
landsatsub4 <- subset(landsat, 5:3)
plotRGB(landsatsub3, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
plotRGB(landsatsub4, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
par(mfrow = c(1, 2))
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
# ?par # help for par() function
landsatsub5 <- subset(landsatcrop, 4:2)
landsatsub6 <- subset(landsatcrop, 5:3)
plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")
# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
plotRGB(landsatsub6, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nFalse Colour Composite")
jpeg(file="landsatRGB-TCC.jpeg")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen
## 2
jpeg(file="landsatRGB-FCC.jpeg")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen
## 2
tiff(file="landsatRGB-TCC.tif")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen
## 2
tiff(file="landsatRGB-FCC.tif")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen
## 2
rf <- writeRaster(landsatcrop, filename = "cropped-landsat.tif", overwrite=TRUE)
rf <- writeRaster(landsatsub5, filename = "cropped-landsat5.tif", overwrite=TRUE)
rf <- writeRaster(s, filename = "cropped-landsat-TCC.tif", overwrite=TRUE)
# S4 method for RasterStackBrick,character
rf <- writeRaster(landsatFCC, filename=file.path("landsatFCC.tif"), format="GTiff", overwrite=TRUE)
rf <- writeRaster(landsatRGB, filename=file.path("landsatTCC.tif"), format="GTiff", overwrite=TRUE)
class(landsatcrop)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"
class(landsat)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
“exploring relationships between raster layers”
## Plot of UV/B reflection “in the ultra-blue wavelength against reflection in the blue wavelength.” “The first plot reveals high correlations between the blue wavelength regions.
Because of the high correlation, we can just use one of the blue bands without losing much information.”
pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")
“This distribution of points between NIR and red is unique due to its triangular shape.
Vegetation reflects very highly in the NIR range than red and creates the upper corner close to NIR (y) axis.
Water absorbs energy from all the bands and occupies the location close to origin.
The furthest corner is created due to highly reflecting surface features like bright soil or concrete.”
pairs(landsatcrop[[4:5]], main = "Red vs NIR")
# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')
# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
## Warning in proj4string(obj): CRS object has comment, which is lost in output
# add the land cover class to the points
ptsamp$class <- over(ptsamp, samp)$class
# extract values with points
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)
## ultra.blue blue green red NIR SWIR1 SWIR2
## [1,] 0.1335017 0.1151333 0.09839138 0.09659141 0.1576821 0.2013586 0.1655977
## [2,] 0.1337186 0.1161743 0.10212145 0.10082026 0.1904286 0.2187511 0.1802577
## [3,] 0.1300753 0.1212923 0.12354765 0.15800740 0.2767190 0.2953259 0.1885419
## [4,] 0.1472076 0.1472509 0.16509888 0.22250289 0.3560046 0.3941944 0.2553578
## [5,] 0.1303572 0.1114032 0.09301314 0.08963006 0.1518918 0.1883684 0.1539304
## [6,] 0.1376221 0.1208368 0.10457201 0.10459370 0.1674193 0.2175801 0.1817541
ms <- aggregate(df, list(ptsamp$class), mean)
# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
## ultra.blue blue green red NIR SWIR1
## built 0.1865284 0.18105759 0.18454494 0.20110170 0.2497486 0.25219162
## cropland 0.1118572 0.08969222 0.08513156 0.05500115 0.4656332 0.15285904
## fallow 0.1325837 0.11699979 0.10439419 0.11388658 0.1784664 0.23141986
## open 0.1392367 0.13802502 0.15345492 0.20746139 0.3433308 0.35600767
## water 0.1340832 0.11714856 0.10011272 0.07951924 0.0493608 0.03392059
## SWIR2
## built 0.20897721
## cropland 0.06948914
## fallow 0.19388931
## open 0.21364008
## water 0.02745309
“The spectral profile shows (dis)similarity in the reflectance of different features on the earth’s surface (or above it).
‘Water’ shows relatively low reflection in all wavelengths, and ‘built’, ‘fallow’ and ‘open’ have relatively high reflectance in the longer wavelengts.”
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
rasterperformed on each pixel / grid cell
raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <-stack(raslist)
landsatRGB <- landsat[[c(4, 3, 2)]]
landsatFCC <- landsat[[c(5, 4, 3)]]
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
to view greenness
Landsat bands NIR = 5, red = 4 #### NDVI one way
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")
vi2 <- function(x, y){
(x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun = vi2)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")
e.g. Identify water
see spectral profile for bands with max / min reflectance
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
Landsat bands: green = 3 and NIR = 5
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Body Index: Landsat Bands Green (3) and NIR (5)", sub = "Monitor water in water bodies")
Landsat bands: NIR = 5 and SWIR = 7
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Index: Landsat Bands NIR (5) and SWIR (7)", sub = "Monitor change in leaf water content")
Landsat bands: NIR = 5 and SWIR = 6 >> not the band to use
ndwi2 <- vi(landsat, 5, 6)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Index: Landsat Bands NIR (5) and SWIR (6)", sub = "Monitor change in leaf water content")
see spectral profile for bands with max / min reflectance
e.g. Identify built areas NDBI (Normalized Difference Built-up Index) https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm Landsat bands: 4 = red, 2 = green, 7 = SWIR
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
ndbi <- bi(landsat, 4, 2, 7)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands red (4), green (2), and SWIR (7)", sub = "Monitor change in built environment (4-2)*7/(4+2)*7")
https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm Landsat bands: 6 = SWIR, 5 = NIR
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
bi <- (bk - bi) / (bk + bi)
return(bi)
}
ndbi <- bi(landsat, 6, 5)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands SWIR (6) and NIR (5)", sub = "Monitor change in built environment (6-5)/(6+5)")
bei <- vi(landsat, 4, 5)
plot(bei, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands NIR (5) and green (3)", sub = "Monitor change in built environment")
par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function
###############
# NDVI plot
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")
###############
# NDWI plot leaf water
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")
###############
# NDWI plot water bodies
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands Green (3) and NIR (5)")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
bsi <- bi(landsat, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")
par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
par(mai = c(0.35, 0.35, 0.35, 0.35))
# ?par # help for par() function
###############
# NDVI plot
ndvi <- vi(landsatcrop, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat NDVI Greeness Index")
###############
# NDWI plot leaf water
ndwi2 <- vi(landsatcrop, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), main = "Leaf Water Index: Landsat Bands NIR (5) and SWIR (7)")
###############
# NDWI plot water bodies
ndwi <- vi(landsatcrop, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), main = "Water Body Index: Landsat Bands Green (3) and NIR (5)")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
bk <- img[[k]]
bi <- img[[i]]
bp <- img[[p]]
bi <- (bk - bi) * bp / (bk + bi) * bp
return(bi)
}
bsi <- bi(landsatcrop, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), main = "Built Index: Landsat Bands Red, Green and SWIR [(4-2)*7/(4+2)*7]")
###############
# Plot Spectra
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")